#include "kernelanalysisimage.h"

kernelAnalysisImage::kernelAnalysisImage()
{

}

kernelAnalysisImage::kernelAnalysisImage(geoLayer* layer,QString savePath)
{
    maxDensity=-9999999;
    minDensity=9999999;
    imageSize=600;
    this->layer=layer;
    this->savePath=savePath;

    getMinMax();
    meanCenter.rx()=(minX+maxX)/2;
    meanCenter.ry()=(minY+maxY)/2;
    //计算格网大小
    double length,lengthX,lengthY;
    lengthX=maxX-minX;
    lengthY=maxY-minY;
    lengthX>lengthY?length=lengthX:length=lengthY;
    gridSize=length/imageSize;
    //初始化格网
    for(int i=0;i<imageSize;i++)
    {
        QVector<double> rows;
        for(int j=0;j<imageSize;j++)
        {
            rows.append(0);
        }
        densities.append(rows);
    }

    foreach (geoObject* object, layer->geoObjects)
    {
        foreach (QPointF coordinate,object->coordinates)
        {
            distances.append(getDistance(meanCenter,coordinate));
        }
    }

    medianDistance=getMedianDistance();
    standardDistance=getStandardDistance();

    double computeDist=medianDistance*(sqrt(1/(log(2))));
    double minDist;
    standardDistance > computeDist ? minDist=computeDist : minDist=standardDistance;
    searchRadius=0.9*minDist*pow(distances.size(),-0.2);
//    searchRadius=2000000;
}

bool kernelAnalysisImage::saveOutPut()
{
    getDensity();
    QImage *image=new QImage(imageSize,imageSize,QImage::Format_RGB32);
    QRgb color = 1;
    for(int row=0;row<imageSize;row++)
    {
        for(int col=0;col<imageSize;col++)
        {
            double normalizedDensity=(densities[row][col]-minDensity)/(maxDensity-minDensity);
            if(normalizedDensity<0.1)   color=qRgb(0,0,0);
            else if(normalizedDensity<0.2)  color=qRgb(80,8,25);
            else if(normalizedDensity<0.3)  color=qRgb(100,10,30);
            else if(normalizedDensity<0.4)  color=qRgb(105,12,38);
            else if(normalizedDensity<0.5)  color=qRgb(110,15,45);
            else if(normalizedDensity<0.6)  color=qRgb(128,19,54);
            else if(normalizedDensity<0.7)  color=qRgb(152,30,60);
            else if(normalizedDensity<0.8)  color=qRgb(199,44,65);
            else if(normalizedDensity<0.9)  color=qRgb(210,50,65);
            else if(normalizedDensity<1)    color=qRgb(238,69,64);
            image->setPixel(row,col,color);
        }
    }

    return image->save(savePath,"jpg",100);
}


double kernelAnalysisImage::getDistance(QPointF pt1, QPointF pt2)
{
    double distance;
    distance=sqrt(pow(pt1.x()-pt2.x(),2)+pow(pt1.y()-pt2.y(),2));
    return distance;
}

void kernelAnalysisImage::getMinMax()
{
    foreach (geoObject* object, layer->geoObjects)
    {
        foreach (QPointF coordinate,object->coordinates)
        {
            if(coordinate.rx()>maxX)    maxX=coordinate.rx();
            if(coordinate.rx()<minX)    minX=coordinate.rx();
            if(coordinate.ry()>maxY)    maxY=coordinate.ry();
            if(coordinate.ry()<minY)    minY=coordinate.ry();
        }
    }
}

double kernelAnalysisImage::getMedianDistance()
{
    double medianDistance;
    std::sort(distances.begin(),distances.end());
    int size=distances.size();
    if(size%2==0)
    {
        medianDistance=(distances[size/2]+distances[size/2+1])/2.0;
    }
    else
    {
        medianDistance=distances[size/2+1];
    }
    return medianDistance;
}

double kernelAnalysisImage::getStandardDistance()
{
    //无权重算法
    double standardDistance;
    double totalX=0,totalY=0;
    foreach(geoObject *object,layer->geoObjects)
    {
        foreach(QPointF coordinate,object->coordinates)
        {
            totalX+=pow(coordinate.x()-meanCenter.x(),2);
            totalY+=pow(coordinate.y()-meanCenter.y(),2);
        }
    }
    int n=layer->geoObjects.size();
    standardDistance=sqrt((totalX+totalY)/n);
    return standardDistance;
}

void kernelAnalysisImage::getDensity()
{
    int size=layer->geoObjects.size();
    for(int row=0;row<imageSize;row++)
    {
        for(int col=0;col<imageSize;col++)
        {
            QPointF pt(gridSize*(row+0.5)+minX,gridSize*(col+0.5)+minY);
            double sum=0;
            for(int i=0;i<size;i++)
            {
                double dist=getDistance(pt,layer->geoObjects[i]->coordinates.first());
                if(dist<searchRadius)
                    sum+=pow(1-pow(dist/searchRadius,2),2);
            }
            double density=(1/pow(searchRadius,2)*3/PI)*sum;
            densities[row][col]=density;
            if(density>maxDensity)  maxDensity=density;
            if(density<minDensity)  minDensity=density;
        }
    }
}
